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Abstract 

o 

' Pressure (density) and velocity boundary conditions inside a flow domain are studied for 2- 

D and 3-D lattice Boltzmann BGK models (LBGK) and new method to specify these conditions 
are proposed. These conditions are consistent with the boundary condition we proposed in Q 
using an idea of bounce-back of non-equilibrium distribution. These conditions give excellent 
results for the regular LBGK models, and were shown to be second-order accurate by numerical 



in 

On 



examples. When they are used together with the improved incompressible LBGK model in 
£jy , the simulation results recover the analytical solution of the plane Poiseuille flow driven by 

I , pressure (density) difference with machine accuracy. 

t: 

O ■ 1 Introduction 

O 

The lattice Boltzmann equation (LBE) method has achieved great success for simulation of trans- 
port phenomena in recent years. Among different LBE methods, the lattice Boltzmann BGK 
model is considered more robust |J. Besides, theoretical discussion is easier for the LBGK due 
to its simple form. Some recent theoretical discussions on LBGK [|], ||, |5| have enhanced our 
understanding of the method and the effect of boundary conditions. In H, analytical solutions 
of distribution functions for plane Poiseuille flow with forcing and plane Couette flow have been 
obtained for the 2-D triangular and square lattice Boltzmann BGK models. It is found that the 
bounce-back boundary condition produces distribution functions with a first-order error compared 
with the analytical distribution functions. In ||, a new technique was developed to seek the ana- 
lytic solution of LBGK model for some simple flow. For example, the velocity profile from the 2-D 
square and triangular LBGK models are shown to satisfy a second-order difference equation of 
the Navier-Stokes equation in the case of plane Poiseuille flow with forcing and Couette flow. The 
technique is generalized in jl]] to include steady-state flows with both x and y velocities, which 
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are assumed to be independent of x. The analysis provides a framework to analyze any velocity 
boundary condition. For example, the analysis explains why the velocity boundary condition for 
the 2-D triangular LBGK model proposed in || generates results of machine accuracy for plane 
Poiseuille flow with forcing. 

In practice, however, a flow is often driven by pressure difference. In general, the pressure 
gradient through the flow field is not a constant and the local pressure gradient is unknown before 
solving the flow. Hence the pressure gradient in many cases cannot be replaced by an external 
force in LBGK computations. In this situation, boundary conditions usually need be implemented 
by giving prescribed pressure or velocity on some "flow boundaries", which are not solid walls 
or interfaces of two distinct fluids. Instead, they are imaginary boundaries inside a flow domain 
(e.g. inlet and outlet in a pipe flow). Their existence is purely for the convenience of study. The 
implementation of these boundary conditions in LBGK is very important but it has not yet been 
well studied. 

Since in lattice Boltzmann method, the pressure is related to the density by the isothermal 
equation of state as p = c 2 s p (c s is the sound speed of the model), a specification of pressure 
difference amounts to a specification of density difference. Early works (see, for example, |7|]) to 
implement pressure (density) flow boundary condition is simply to assign the equilibrium distri- 
bution computed with the specified density and some velocity (maybe zero) to the distribution 
function. This method introduces significant errors: the real pressure gradient obtained in the 
simulation for the Poiseuille flow is not a constant. It is approximately a constant only some dis- 
tance away from the inlet and outlet of the channel. Besides, even away from the inlet and outlet 
region, the pressure gradient is different from the intended value. Maier et al. ||] proposed an 
alternative pressure or velocity flow boundary condition for the 3-D 15-velocity direction LBGK 
model, and their results are greatly improved over the equilibrium distribution approach. The 
pressure or velocity flow boundary condition in is obtained through a post-streaming rule to 
the distribution functions based on an extrapolation. However, this pressure or velocity boundary 
condition is still to be improved due to some inconsistency (see discussion in Section 2). Its inac- 
curacy is most noticeable in the following case: when this pressure or velocity boundary condition 
is applied to the modified LBGK [|[, which corresponds to a macroscopic momentum equation 
having the analytical solution of Poiseuille flow with pressure (density) gradient, the simulation 
results are far from the analytical results. 

In this paper, we propose a general way to specify pressure or velocity on flow boundaries. The 
implementation is a natural extension of the wall boundary condition described in our previous 
paper Q . The result shows a clear improvement to the flow boundary conditions in |8| for ordinary 
LBGK models. Besides, for the modified LBGK model, These flow boundary conditions produce 
results of machine accuracy for Poiseuille flow. 

2 Pressure or Velocity Flow Boundary Condition of the 2-D 
Square Lattice LBGK Model 

2.1 Governing Equation 

The square lattice LBGK model (d2q9) is expressed as (H,fli~OH,|llll): 

f i ( x + 5e i ,t + 5)-f t (x,t) = -kfi(x,t)-ft 9 \x,t)}, i = 0,l,...,8, (1) 
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where the equation is written in physical units. Both the time step and the lattice spacing have 
the value of 8 in physical units. /j(x, t) is the density distribution function along the direction ej 
at (x,i). The particle speed e^'s are given by ej = (cos(-7r(i — l)/2), sin(7r(z — l)/2),i = 1,2,3,4, 
and ej = \/2(cos(7r(z — 4 — i)/2), sin(-7r(i — 4 — i)/2),z = 5,6,7,8. Rest particles of type with 
eg = is also allowed (see Fig. 1). The right hand side represents the collision term and r is the 
single relaxation time which controls the rate of approach to equilibrium. The density per node, 
p, and the macroscopic flow velocity, u = (u x ,u y ), are defined in terms of the particle distribution 
function by 

8 8 

fi = p> J2 fo ei = p u - ( 2 ) 

i=0 i=l 

The equilibrium distribution functions f^ eq \x, t) depend only on local density and velocity and 
they can be chosen in the following form (the model d2q9 |To|j): 



ft q) = ^[l + 3(e J .u) + ^(e,.u) 2 -^u-u], i = 1,2,3,4 (3) 

ft q) = ^p[l + 3(e,.u) + ^( ei .u) 2 -^u.u], z = 5,6,7,8. 

A Chapman-Enskog procedure can be applied to Eq. (JxJ) to derive the macroscopic equations of 
the model. They are given by: the continuity equation (with an error term 0(5 2 ) being omitted): 

^ + V-(pu)=0, (4) 

and the momentum equation (with terms of 0{5 2 ) and 0(5u 3 ) being omitted): 

dt(pu a ) + dpipuaup) = -d a (c 2 s p) + df3(2vpS al3 ), (5) 

where the Einstein summation convention is used. S a p = h(d a up + dpu a ) is the strain-rate tensor. 

9 2 1 2r - 1 

The pressure is given by p = ctp, where c s is the speed of sound with c~ = — , and v = o, 

3 6 

with v being the the kinematic viscosity. The form of the error terms and the derivation of these 



equations can be found in [^, 13]. 



In this paper, we will take the Poiseuille flow as an example to study the pressure (density) or 
velocity inlet/outlet condition. The analytical solution of Poiseuille flow in a channel with width 
2L for the Navier-Stokes equation is given by: 

y 2 dp dp 

u x = u {l-j^), u y = 0, ^ = -G, ^ = °> ( 6 ) 

where the pressure gradient G is a constant related to the centerline velocity uq by 

G = 2puu /L 2 , (7) 

and the flow density p is a constant. The Reynolds number is defined as Re = Uq(2L)/v. 

The Poiseuille flow is an exact solution of the steady-state incompressible Navier-Stokes equa- 
tions: 

V • u = 0. (8) 
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p 

dp(ll a Up) = -da(-) + vdp(]U a , (9) 

Po 

On the other hand, the steady-state macroscopic equations of LBGK model, Eq. ([!]), is given by: 

V • (pu) = 0, (10) 
2r — 1 

dp(pu a Uf3) = -d a (c 2 s p) H - — 5dpdp(pu a ). (11) 

These equations are different from the incompressible Navier-Stokes equations Eqs. (||J9|) by terms 
containing the spatial derivative of p. These discrepancies are called compressibility error in LBE 
model. Thus, the Poiseuille flow given by Eq. (|6|) is not the exact solution of Eqs. ( |10| , |TI| ) when 
pressure (density) gradient drives the flow. That is, due to change of pressure (density) in the 
x— direction, u x is not constant in the x— direction, and the velocity profile of the solution of 
Eqs. fli~0| , 11) is no loger a parabolic profile. For a fixed Mach number (mq fixed), as 5 —* 0, 



the velocity of the LBGK simulation will not converge to the velocity in Eq. (q) because the 
compressibility error becomes dominating. This phenomenon is seen in the result in (8|, where 
the error of velocity increases as the number of lattice grid increases. Besides, from d x {pu x ) = 
(suppose u y = in the simulation), one can see that u x should be increasing linearly along the 
flow direction since p is decreasing linearly. This makes the comparison of u x with the analytical 
velocity of Poiseuille flow somehow ambiguous. 

To make a more sensible study for Poiseuille flow with pressure (density) or velocity flow 
boundary condition, it is better to use the improved incompressible LBGK model proposed in 
[||]. The model (called d2q9i) is given by Eq. @) with the same ej and the following equilibrium 
distributions: 

Am) 4 r 3 i 

ft q) = ^[p + 3e i .v + ^(e l .v) 2 -^vv], i = l,2,3,4, 

ft Q) = ^ + 3e ! .v+^(e l .v) 2 -^.v] )l = 5,6,7,8, (12) 



and 



E /< = E ft Q) = p, E fa = E f^^ = (13) 



where v = (v x ,v y ) (like the momentum in the ordinary LBGK model) is used to represent the 
flow velocity. The macroscopic equations of d2q9i in the steady-state case (apart from error terms 
ofO(5 2 ): 

V-v = 0, (14) 

d/3{v a vp) = -d a (c 2 s p) + vdppVa, (15) 

are exactly the steady-state incompressible Navier-Stokes equation. In the model, pressure is 
related to density by (? s p = p/po {c 2 s = 1/3) for a flow with constant density like Poiseuille flow, 
and v = 2r r 1 ^. 

If the flow is steady, d2q9i is superior to d2q9 to simulating incompressible flows. If the flow 
is unsteady, one may consider the continuity equation derived from d2q9 given by dfp+ (V • p)u+ 
pV • u = 0. In a situation where the first two terms likely cancel each other, d2q9 may be of 
advantage (to approximate the continuity equation V • u). If the first two terms have the same 
sign, then d2q9i is better. 
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2.2 Review of The Velocity Wall Boundary Condition 

It is proved in Q || that if the flow is steady and independent of x, then the solution /, of 
Eq. (|H) produces a velocity profile that satisfies a difference equation which is a second-order 
approximation of the Navier-Stokes equation. If the boundary condition is chosen correctly, then 
the difference equation near the boundary is consistent with the difference equation inside. 

A velocity wall boundary condition is proposed in ||] as follows: take the case of a bottom 
node in Fig. 1, the boundary is aligned with x— direction with f&, h, /§ pointing into the wall. 
After streaming, /o, fi, hi Ja, hi h are known. Suppose that specified on the wall, we 

need to determine hi hi h and P from Eqs. (|2j), which can be put into the form: 

h + h + k = P - (/o + fi + h + h + h + h), (16) 

h-h = pu x -(h-h-h + h), (17) 
h + h + h = pu y + (h + h + h)- (is) 

Consistency of Eqs. (16|,1^) gives 

P=—^—lfo + fi + h + 2(h + h + h)}. (19) 
1 — u,, 



We assume the bounce-back rule is still correct for the non-equilibrium part of the particle 
tribution norr 
can be found as 



distribution normal to the boundary (in this case, h ~ = h ~ fi)- With h known, hi h 



H = U + \tpUyi 

h = h ~^(fi- h) + ^pux + ^pu y , 

h = h + ^(fi- h) -^pux + ^puy. (20) 

The collision step is applied to the boundary nodes also. For non-slip boundaries, this boundary 
condition is reduced to that in [Hi. 



2.3 Specification of Pressure on a Flow Boundary 

Now let us turn to pressure (density) flow boundary condition. Its derivation is based on Eq. (||) 
as for velocity wall boundary condition. Suppose a flow boundary (take the inlet in Fig. 1 as 
example) is along the y— direction, and the pressure is to be specified on it. Suppose that u y is 
also specified (e.g. u y = at the inlet in a channel flow). After streaming, hi h, f4i hi h are 
known, p = pi n , u y = are specified at inlet. We need to determine u x and fi, hi h from Eq. (||) 
as following: 

h + h + h = Pin - (h + h + h + h + h + h), (21) 

h + h + h = p in u x + (/ 3 + h + h) i (22) 

h- h = h- h + h- h- (23) 
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Consistency of Eqs. ( ^Hf2^ ) gives 



1 _ [/0 + /2 + /4 + 2(/ 3 + /6 + /7)] (M) 



We use bounce-back rule for the non-equilibrium part of the particle distribution normal to 
: inlet tc 
equations: 



the inlet to find fa — /{ = fa — '. With fa known, fa, fa are obtained by the remaining two 



fa = fa + 7^PinU x , 

fa = fa ~\{f2- fl) + \pinU x , 

fa = k + \(f2~ h) + \p in u x . (25) 

The corner node at inlet needs some special treatment. Take the bottom node at inlet as 
an example, after streaming, fa, f 4, fa are known; p is specified, and u x = u y = 0. We need to 
determine fa, fa, fs, fa, fs- We use bounce-back rule for the non-equilibrium part of the particle 
distribution normal to the inlet and the boundary to find: 

fi = h + {f[ eq) - ft 9) ) = fa, fa = h + {f[ eq) - fi eq) ) = h, (26) 

Using these /i, fa in Eqs. (||), we find: 

h = h, k = h = \\Pin - (/o + h + h + h + U + h + fi)\- (27) 

Similar procedure can be applied to top inlet node and outlet nodes including outlet corner nodes. 
The case that p and non-zero u y is specified at a flow boundary along y— direction can be handled 
in the same way. 

Here, it is useful to compare our pressure boundary condition with that proposed in []|], which 
is given by the following post-streaming rule (an extrapolation) at inlet: after streaming, fi, fs, fs 
are calculated as 

/ f (x,t) = /+(x,t-5)-(/ i (x,t)-/+(x,t-,5)), z = 1,5,8 (28) 

where /j + (x, t — 5) = /j(x, t — 5) — ^(/j(x, t — 5) — (x, t — 5)) is the distribution functions at 
previous time step after collision and before streaming, fj is along ej with = — 2e n (the 
inner normal e n = e± in the case). Thus, for i = 1,5, 8, j = 3, 6, 7 respectively. The density p is 
set to the specified inlet value and u y is set to zero to compute f\ eq \~x.,t). At the bottom, f\, fs 



are computed using Eq. ( |28| ) and then f2,fs,fe are obtained in the treatment of wall boundary 
condition. Notice, however, that at the inlet, Ya=o fi ma y n °t be equal to the specified density 
and Y^i=i e iyfi ma Y n °t De equal to zero with this post-streaming operation. This inconsistency 
causes some inaccuracy in simulations and leaves room for improvement. 

2.4 Specification of Velocity on a Flow Boundary 

In some calculations, velocities u x ,u y are specified at a flow boundary (take the inlet in Fig. 1 as 
example). In the case of channel flow, after streaming, fz, f&, fz, f§, fa are known at inlet. u x ,u y 
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are specified at inlet (for the special case of Poiseuille flow, u y = 0), we need to determine p and 
fa, f 5, fa- This is actually equivalent to a velocity wall boundary condition. Using our velocity 
wall boundary condition in previously described, we find: 

P = T ^—[fa + fa + fa + 2(fa + fa + fa)), 
l-u x 

fa = fa + ^pu x , 

fa = fa ~ 2C/2 - fa) + \p u y + qPUx, 

h = fa + 7^/2 - fa) ~ ^pu y + ^pu x . (29) 

The effect of specifying velocity at inlet is similar to specifying pressure (density) at inlet. Density 
difference in the flow can be generated by the velocity inlet condition. 

At the inlet bottom (non-slip boundary), special treatment is needed. After streaming, 
/l) fa, fa} fa, fa need to be determined. Using bounce-back on normal distributions gives: 

/1 = fa, fa = fa, 

expressions of x, y momenta give: 

fa - fa + fa = -(fa -fa- fa) = fa, 

fa + fa-fa = -(fa-fa-fa) = fa, (30) 

which gives 

fa = fa, 

fa = fa = \[p~ (fa + fa + fa + fa + fa + fa + fa)}, (31) 

but there is no more equation available to determine p. The situation is similar to a corner wall 
node (the intersection of two perpendicular walls). In the situation, p at the inlet bottom node 
can be taken as the p of its neighboring flow node, thus the velocity inlet condition is specified. 

From the discussion given above, we can unify boundary conditions (on a wall boundary or in 
a flow boundary) in 2-D simulation on a straight boundary as: 

• Given u x ,u y , find p and unknown /j's. 

• Given p and the velocity along the boundary, find the velocity normal to the boundary and 
unknown /j's. 

Eq. (||) is used to determine p or the normal velocity and the unknown /j's. The formula are 
given in sections 2.2, 2.3, 2.4. 

Again, it is useful to compare our velocity flow boundary condition with that proposed in ||, 
which is given by the following post-streaming rule (a zeroth-order extrapolation) at inlet: after 
streaming, fa, fa, fa are calculated using 

/ f (x,t) = /+(x,t-5), z = 1,5,8 (32) 
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where f^(x.,t — 5) is the distribution function at previous time step after collision. After p is 

computed using /j's, f\ eq \x.,t) can be computed using this p and the specified velocities u x ,u y . 
In this approach, the determination of unknown /j's does not use the information of known /j's at 
present time. This is inconsistent with the present distribution in the flow. Suppose that initially, 
fl eq \ i = 0, • • • , 8 are computed by using some density po and the specified u x ,u y , and one assigns 
fi = f\ eq \ i = 0, • • • ,8, then collision does not change fi, and the post-streaming rule Eq. (|3 

does not change fi and p. Hence fi = // for all time in the simulation. This velocity inlet 
condition amounts to assign the equilibrium distribution to fi and it makes a significant error. 
The result is worse than that of pressure inlet condition ||. 

2.5 Boundary Conditions for the Modified Incompressible Model d2q9i 

The velocity wall boundary condition and flow boundary conditions for d2q9i are similar to that 
of d2q9. It is from equations J2i=o fi = P an d J2i=i e ifi = v an d hence some modifications are 
needed as follows: 



In wall boundary condition, Eq. (|19| ) is replaced by 

P = v y -[f + h + h + 2(/ 4 + f 7 + h)]. (33) 

and in Eq. (p0|), pu x ,pu y are replaced by v x ,v y respectively. 



In pressure flow boundary condition, Eq. (24) is replaced by 



v x = p- [f + h + h + 2(/ 3 + h + h)\ , (34) 



and in Eq. (^5[), pi n u x is replaced by v x . 



Similar replacement in velocity flow boundary condition, Eq. (|29|). 



3 Numerical Results and Discussion 

We report and discuss the numerical results for Poiseuille flow with pressure (density) or velocity 
flow boundary condition. The simulation is performed on both models d2q9 and d2q9i. The main 
result in the simulation of d2q9i is the achievement of machine accuracy. The main result in 
the simulation of d2q9 is the achievement of second-order accuracy of the boundary conditions. 
The width of the channel is assumed to be 2L = 2. We use nx, ny lattice nodes on the x— and 
y— directions, thus, 5 = 2/{ny — 1). The initial condition is to assign fi = f^ computed using a 
constant density po, and zero velocities. The steady-state is reached if 



EiEj \u x {i,j,t + 5) - u x (i,j,t)\ + \u y (i,j,t + 5) - u y {i,j,i)\ 



< 5 • Tol. 



EiEj \ux(i,j,t)\ + \u y (i,j,t)\ 

For model d2q9i, u x ,u y are replaced by v x ,v y . Tol is a tolerance usually set to 10 -10 . On the 
wall, boundary condition discussed in section 2.2 is used to make non-slip boundaries. 
We also define a LI error as: 

_ £t£jl*4foj) -Ux(i,j)\ + \ul(i,j) -u y (i,j)\ 
where u x ,u y is the analytical velocity. 
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3.1 Results of Model d2q9i 



For model d2q9i, we carried out simulations with a variety of Re, nx, ny, uq (uq is the peak velocity 
in the channel) using the pressure or velocity flow boundary condition. The range of Re is from 
0.0001 to 30.0; the range of r is from 0.56 to 20.0 and the range of uq is from 0.001 to 0.4; the 
largest density difference simulated (not the limit) is pi n = 5.6, p ou t = 4.4 with nx = 5,ny = 3 
corresponding to a pressure gradient of G = 0.1. The magnitude of average density po is irrelevant 
for the simulation Q. 

For all cases where the simulation is stable, the steady-state velocity and density show: 

• The velocity field v x is accurate up to machine accuracy compared to the analytical solution 
in Eq. (^), v y is very small with maximum of \v y \ in the whole region being in the order of 
10~ 12 . For example, for nx = 5,ny = 3, uq = 0.1, r = 0.56, Re =10, the relative LI error of 
v in the whole flow region is 0.485 • 10~ 10 . 

• The density is uniform in the cross channel direction, and linear in the flow direction. The 
density difference p(i + — p(i,j) is a constant through the flow region, its value is equal 
(up to machine accuracy) to the analytical value set by the constant pressure gradient. 

• Velocity v x is uniform in the x— direction, the results are the same for different nx. 

If the computed velocity were plotted with the analytical velocity, there would be no difference 
to naked-eyes. 

It is also noticed that with pressure (density) gradient to drive Poiseuille flow, the maximum 
Reynolds number which makes the simulation stable is far less than that with external forcing. 
For ny = 5, the maximum Re is about 30. Refinement of mesh can increase the maximum Re. 

When the pressure flow boundary condition is replaced by the method in (8|, machine accuracy 
can no longer be obtained, for a simple case nx = 17, ny = 9, uo = 0.03542, r = 0.67, Re =5, with 
a moderate pressure gradient of G = 0.001004, pi n = 5.006, p ou t = 4.994, the relative LI error of 
v in the whole flow region is 0.1824 • 10 -2 , with maximum of \v y \ being 0.1364 • 10~ 3 , not very 
small compared to uq. The density difference (p(i + 1, j) — p(i,j))/5 is no longer a constant, its 
range is from -0.002094 to -0.003872 (the analytical value is -0.003012 = -G/c 2 s = -3G). The 
result indicates that the pressure or velocity flow boundary condition proposed in this paper is a 
clear improvement of that in ||. 

Similar results of machine accuracy are obtained by specifying the analytical velocity profile 
given in Eq. ([]) at inlet and pressure (density) at outlet by using the flow boundary conditions in 
this paper. In the case, there is a uniform pressure (density) difference in the region. The value 
of the density difference depends on no and the outlet density. 

3.2 Results of Model d2q9 

Since the ordinary LBGK model d2q9 is still widely used for simulations. It is worthwhile to 
do some simulations with d2q9 with our flow boundary conditions and show that they give a 
second-order accuracy. We use d2q9 to Poiseuille flow with pressure or velocity flow boundary 
condition. Since as 5 — > the computed solution does not approach the analytical solution, we 
will use the result of the finest mesh as the exact solution to compute the error. Simulations with 
successively doubled lattice steps are carried out to observe the convergence. The example uses 
fixed Re = 10, uq = 0.1, po = 5. The pressure gradient G in Eq. @ and then pressure (density) 
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at inlet /outlet can be obtained as G = 0.02 and pi n = 5.12, p out = 4.88 respectively to be used 
in the pressure (density) flow boundary condition. For the velocity flow boundary condition, the 
analytical velocities in Eq. @ are used at inlet, and po is specified at outlet. Thus, the results 
of the pressure flow boundary condition and the velocity flow boundary condition are similar but 
not identical. Define Ix = nx — 1, ly = ny — 1 to represent the number of lattice steps in x— 
and y— directions, we use Ix = 4, 8, 16, 32, 64, 128, 256, ly = 2, 4, 8, 16, 32, 64, 128 respectively to do 
the simulation. The LI error is defined in Eq. ( |35[ ) with u x ,u y being replaced by the computed 
velocities of Ix = 256, ly = 128. r has to be changed as Ix, ly are changed to keep the same Re, no 
(values of r are included in Table I). The convergence result is summerized in Table I. The case 
of Ix = 4, ly = 2 was not shown in Table I, because the simulation is unstable for the velocity 
inlet condition. The ratio of two consecutive LI errors is also shown. The ratio is approximately 
equal to 4, indicating a second-order accuracy. For all runs, it is also observed that 

• Velocity u x is monotonically increasing in the x— direction, the result is not sensitive to 
the value of Ix, we usually take Ix = 2 ly. If the pressure boundary condition in ^ is 
used, u x on the center line of the channel decreases first, then increases, and decreases again 
near the outlet. This behavior deviates from the macroscopic continuity equation of LBGK 
d x {pu x ) = in the case, indicating that errors are introduced by the pressure boundary 
condition in M. An example of centerline velocity function of x is presented in 
Fig. 4. 

• u y is very small compared to no, with maximum of \u y \ being approximately of 10 -6 • no in 
typical cases. If the pressure boundary condition in Q is used, maximum value of \u y \ is 
like 10~ 3 • no typically. 

• The density is uniform in the cross channel direction, and linear in the flow direction. The 
density difference p(i + 1, j) — p(i, j) is almost a constant through the flow region, its value is 
equal to the analytical value with a fluctuation of less than 0.03 % in a worst case observed 
with Re= 5. If the pressure boundary condition in || is used, The fluctuation of density 
difference p(i + 1, j) — p(i,j) from the analytical value reaches 20 % for the same case 
mentioned above with Re = 5. 

In the case where the density gradient is small, the computed velocity profile are close to the 
analytical velocity profile of Poiseuille flow. We present an example here with nx = 9, ny = 5, 
Re=0.8333, n = 0.008333 , G = 0.001667, p in = 5.01, p out = 4.99 with both pressure (density) 
flow boundary conditions in this paper and in M. Fig. 2 shows the velocity u x as a function of y 
at i = 5 (the middle section of the channel), Fig. 3 shows the centerline density profile along the 
x direction, and Fig. 4 shows the centerline u x along the x direction, the solid line represents the 
corresponding analytical solution and the symbols o, + represent the computed solutions with the 
pressure boundary conditions in this paper and in Q respectively. Both computed velocities u x 
at the mid-channel has no difference to the analytical solution to naked eyes (the relative error 
are of order 10 -3 for both pressure boundary conditions). The computed centerline density with 
our method looks identical to the analytical solution (given as a linear function crossing pi n and 
p ou t at inlet and outlet respectively), while the centerline density with the method in ||] has a 
discernible difference with the linear function especially near the inlet and outlet. The centerline 
u x with our method is monotonically increasing, the behavior is consistent with the continuity 
equation of LBGK d x (pu x ) = in the case, while the centerline u x with the method in f|] has a 
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behavior inconsistent with the continuity equation near the inlet and outlet. This again shows an 
clear improvement of our method to that in Q. 



4 Flow Boundary Conditions and Preliminary Results for the 
3-D 15-velocity LBGK Model 

Since 3-D model is needed in practical problems, we give a discussion of the pressure or veloc- 
ity flow boundary condition for the 3-D 15-velocity LBGK model (d3ql5) and present a brief 
statement about its simulation results. The model is based on the LBGK equation Eq. (|l]) with 
£ = 0, 1, • • • , 14, where e$, £ = 0, 1, • • • , 14 are the column vectors of the following matrix: 



E 



1-10 
1 





-1 
1 



1 
1 

-1 1 



-1 1 -1 
1 -1 1 
-1 -1 1 



and ej, £ = 1, ■ ■ ■ , 6 are clasified as type I, e$, i = 7, ■ • ■ , 14 are clasified as type II. The density per 
node, p, and the macroscopic flow velocity, u — {^xi 

u y ,u z ), are defined in terms of the particle 

distribution function by 

14 14 

p, fc ei = p u - 



1=0 



(36) 



i=l 



The equilibrium can be chosen as: 



e( e <l) 



(eq) 



1 1 

- 8 P ~ ~ 3 PU • u, 

-p + -pe, ■ u + -pie, ■ u) -- 



(eg) 



1 



1 



1 



pu • u, 

1 



777 P + 777 P e i " u + 777/°( e i • u) - — pu • U. £ = II 



64 



24 



16 



48' 



(37) 



Since we will use this model for 2-D simulation in this paper, it is clear to give a projection 
of the velocities in the xz plane as shown in Fig. 1. The y— axis is pointing into the paper, so are 
velocity directions 3,7,9,12,14, while the velocity directions 4,8,10,11,13 are pointing out (velocity 
directions 3,4 have a projection at the center and are not shown in the figure). The flow direction 
is still x, and the cross channel direction is z. The macroscopic equations of the model is the 
same as Eqs. (|,|) with c 2 s = 3/8, and v = (2r - 1)5/6. 

The velocity wall boundary condition proposed in |j|] has the following version for the model 
d3ql5: take the case of a bottom node (wall node) as shown in Fig. 1, the wall is on xy plane. 
After streaming, /j, (£ = 0, 1, 2, 3, 4, 6, 8, 9, 12, 13) are known, Suppose that specified 
on the wall, we need to determine fi,i = 5,7,10,11,14 and p from Eqs. (j36j). Similar to the 
derivation in d2q9, p is determined by a consistency condition as: 



1 



The expression of z— momentum gives: 



[/o + fi + h + h + h + 2(/ 6 + h + h + fi2 + /: 



13 J 



(38) 



h + fi + flO + hi + /l4 = PU Z + (h + h + h + fl2 + /l3), 



(39) 
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If we use bounce-back rule for the non-equilibrium part of the particle distribution fi,(i = 
5,7, 10, 11, 14) to set 

fi = fi + l + (ti eq) -fk g h i = 5,7,ll 

fi = fi-i + (ri eq) -f£ q h i = 10,14 (40) 

then Eq. (j39|) is satisfied, and all fi are defined. In order to get the correct x— , y— momenta, 
we further fix this f$ (bounce-back of non-equilibrium fi in the normal direction) and modify 
/7,/io,/n,/i4 as in §]: 

11 

fi<-fi + je ix 6 x + -e iy S y . i = 7, 10, 11, 14 (41) 

This modification leaves z— momentum unchanged but adds 5 x ,5 y to the x—,y— momenta re- 
spectively. A suitable choice of 5 X and 5 y then gives the correct x— , y— momenta. Finally, we 
find: 

h = h + T^puz, 

l l 
12 



fi = fj + t^P u z + -j[ e ix(pu x - fi + h) + e iy (pu y - h + / 4 )], (42) 



where j is the index corresponding to = — ej (e.g., j = 8 for i = 7 and j = 9 for i = 10). In 
the case of non-slip boundary, this boundary condition is reduced to that in Q. 

The derivation of pressure (density) flow boundary condition uses a similar way as for velocity 
wall boundary condition. Suppose the boundary (take the case of inlet in Fig. 1) is on yz plane with 
specified pi n and u y = u z = 0. After streaming, we need to determine u x and fi, (i = 1, 7, 9, 11, 13) 
from Eq. (36). The consistency condition gives: 

PinU x — Pin [/o + h + h + h + h + 2(/ 2 + h + fw + fn + /u)] , (43) 

which determines u x at inlet, using a similar procedure as in deriving the boundary condition, we 
find: 

fl = h+ 7^PinU x , 

fi = fj + jxPinUx ~ ^hy(f3 ~ fi) +e iz (f 5 - /e)], i = 7,9,11,13 (44) 

where j direction is opposite to i direction. 

Same procedure as in d2q9 is used at the inlet bottom node (non-slip) to derive: 

fl = /2, h = /6) 

1 1 

fl = f&—-Z\h — U)i hi = /l2 + — A)) 

/9 = /l0 = /l3 = /l4 

= \[pxn ~ (fo + h + h + h + Ia + k + h + h + h + fn + fn)}, (45) 
and similar results for inlet top and outlet condition. 
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The pressure (density) boundary condition in [|| is specified as post-streaming rule (take the 
inlet as an example): 

fi(x,t) = f?(x,t-6)-(f j (x,t)-f+(xi,t-S)), { = 1,7,9,11,13 (46) 

where f^(x,t — 5) is the distribution functions at previous time step after collision, fj is along 
ej with ej = — 2e n (the inner normal e n = ei in the case). Thus, j = 2,14,12,10,8 for 
i = 1,7,9,11,13 respectively. Then p is set to pi n and u y ,u z are set to zero to compute . 
Again, the density and z— momentum from fi may not be correct, indicating some inconsistency. 

The velocity flow boundary condition, which can be viewed simply a velocity wall boundary 
condition, can be derived similarly. But the corner node with non-slip condition needs some 
treatment as in the d2q9 case. The details are easy to work out and omitted here. 

Simple modifications are used to derive wall boundary condition, pressure (density) or velocity 
flow boundary condition for the improved incompressible model d3ql5i, which is the counterpart 
of d2q9i in 3-D case. 

Simulations on plane Poiseuille flow are performed on d3ql5 and d3ql5i using the pressure or 
velocity flow boundary condition. The only difference with 2-D simulations is periodic condition 
is used on y— direction, and initial condition is uniform in y direction with zero velocity. It is 
found that the result is uniform in the y-direction and is independent of the number of nodes in 
y— direction. u y is identically zero (in the order of 10 -16 because of round-off error). The results 
are very similar, although not identical, to the results of d2q9, d2q9i. Again, d3ql5i with our 
boundary condition and pressure (density) flow boundary condition gives results with machine 
accuracy, showing a clear improvement over the pressure boundary condition in [p. 
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Table I. LI relative error of velocities as lx,ly are doubled, Re=10, uo = 0.1 with our pres- 
sure or velocity flow boundary condition. Error are from comparison with the velocities of 
Ix = 256, ly = 128 (r = 4.34). The symbol (-2) represents 10 -2 . Ratio of two consecutive 
LI errors is also shown. 





lx 


8 


16 


32 


64 


128 


inlet 


ly 


4 


8 


16 


32 


64 


condition 


r 


0.62 


0.74 


0.98 


1.46 


2.42 


pressure 


LI error 


0.1049(-2) 


0.2522(-3) 


0.6135(-4) 


0.1458(-4) 


0.2915(-5) 


(density) 


ratio 


4.159 


4.110 


4.208 


5.003 




velocity 


LI error 


0.2301(-3) 


0.4882(-4) 


0.1167(-4) 


0.2774(-5) 


0.5582(-6) 




ratio 


4.713 


4.183 


4.207 


4.970 
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6 Figure Caption 

Fig. 1, Schematic plot of velocity directions of the 2-D (d2q9) model and projection of 3-D 
(d3ql5) model in a channel. In the 3-D model, The y— axis is pointing into the paper, so are 
velocity directions 3,7,9,12,14, while the velocity directions 4,8,10,11,13 are pointing out (velocity 
directions 3,4 have a projection at the center and are not shown in the figure). 



Fig. 2, Velocity function of y at i = 5 (the middle section of the channel) for the 

case nx = 9,ny = 5, Re=0.8333, u = 0.008333 , G = 0.001667, p in = 5.01, p out = 4.99 with 
pressure (density) flow boundary conditions. The solid line represents the corresponding analytical 
solution of Poiseuille flow. The symbols o, + represent the computed solutions with the pressure 
flow boundary conditions in this paper and in H respectively. 



Fig. 3, Centerline density profile along the x direction for the case nx = 9, ny = 5, Re=0.8333, 
uq = 0.008333 , G = 0.001667, pi n = 5.01, p ou t = 4.99 with pressure (density) flow boundary con- 
dition. The solid line represents the corresponding analytical solution (a linear function crossing 
Pin and p ou t at inlet and outlet respectively). The symbols o, + represent the computed solutions 
with the pressure flow boundary conditions in this paper and in Q respectively. 



Fig. 4, Centerline u x along the x direction for the case nx = 9, ny = 5, Re=0.8333, no = 
0.008333 , G = 0.001667, pi n = 5.01, p m t = 4.99 with pressure (density) flow boundary condition. 
The solid line represent the analytical solution of Poiseuille flow. The symbols o, + represent the 
computed solutions with the pressure flow boundary conditions in this paper and in || respec- 
tively. 
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